## results_appendixfigures.R
## Vito D'Orazio and Idean Salehyan
## replication of figures in the appendix comparing base models to those with full set of controls and those with entropy balancing
## use this script to replicate: figure A1, figure A2, figure A3, figure A4, figure A5, figure A6

rm(list=ls())

library(foreign)
library(readstata13)
library(ggplot2)
library(gridExtra)
library(MASS)

# to replicate results, set this to your working directory
setwd("/Users/vjdorazio/Desktop/TerrorExperiment/paper/II/Final/replication")

mydata <- read.dta13("analysisdata2.dta")
mydata$webal <- mydata$"_webal"

mydata2 <- read.dta("analysisdataJune20.dta")


## replicate model 1
fit.1 <- glm(terror ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")

## model 1 with control dummies
fit.1controls <- glm(terror ~ group + white + arab + white_group + arab_group + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata, family="binomial")
out.1controls <- summary(fit.1controls)

## model 1 with entropy balancing
fit.1eb <- glm(terror ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial", weights=webal) # warning is ok
out.1eb <- summary(fit.1eb)

## replicate model 2
fit.2 <- glm(shooting ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")

## model 2 with control dummies
fit.2controls <- glm(shooting ~ group + white + arab + white_group + arab_group + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata, family="binomial")
out.2controls <- summary(fit.2controls)

## model 2 with entropy balancing
fit.2eb <- glm(shooting ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial", weights=webal) # warning is ok
out.2eb <- summary(fit.2eb)

#### Table 3 ####

## replicate model 3
fit.3 <- glm(shooting ~ control + white , data=mydata[which(mydata$all_group==0),], family="binomial")

fit.3controls <- glm(shooting ~ control + white + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata[which(mydata$all_group==0),], family="binomial")
out.3controls <- summary(fit.3controls)

fit.3eb <- glm(shooting ~ control + white , data=mydata[which(mydata$all_group==0),], family="binomial", weights=webal)
out.3eb <- summary(fit.3eb)

## replicate model 4
fit.4 <- glm(terror ~ control + white, data=mydata[which(mydata$all_group==0),], family="binomial")

fit.4controls <- glm(terror ~ control + white + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata[which(mydata$all_group==0),], family="binomial")
out.4controls <- summary(fit.4controls)

fit.4eb <- glm(terror ~ control + white, data=mydata[which(mydata$all_group==0),], family="binomial", weights=webal)
out.4eb <- summary(fit.4eb)

## replicate model 5
fit.5 <- glm(shooting ~ group + white_group, data=mydata[which(mydata$all_group==1),], family="binomial")

fit.5controls <- glm(shooting ~ group + white_group + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata[which(mydata$all_group==1),], family="binomial")
out.5controls <- summary(fit.5controls)

fit.5eb <- glm(shooting ~ group + white_group, data=mydata[which(mydata$all_group==1),], family="binomial", weights=webal)
out.5eb <- summary(fit.5eb)

## replicate model 6
fit.6 <- glm(terror ~ group + white_group, data=mydata[which(mydata$all_group==1),], family="binomial")

fit.6controls <- glm(terror ~ group + white_group + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata[which(mydata$all_group==1),], family="binomial")
out.6controls <- summary(fit.6controls)

fit.6eb <- glm(terror ~ group + white_group, data=mydata[which(mydata$all_group==1),], family="binomial", weights=webal)
out.6eb <- summary(fit.6eb)

#### Table 4 ####

## replicate model 7
fit.7 <- glm(motive_criminal ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")

fit.7controls <- glm(motive_criminal ~ group + white + arab + white_group + arab_group + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata, family="binomial")
out.7controls <- summary(fit.7controls)

fit.7eb <- glm(motive_criminal ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial", weights=webal)
out.7eb <- summary(fit.7eb)

## replicate model 8
fit.8 <- glm(motive_racial ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")

fit.8controls <- glm(motive_racial ~ group + white + arab + white_group + arab_group + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata, family="binomial")
out.8controls <- summary(fit.8controls)

fit.8eb <- glm(motive_racial ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial", weights=webal)
out.8eb <- summary(fit.8eb)

## replicate model 9
fit.9 <- glm(motive_religious ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")

fit.9controls <- glm(motive_religious ~ group + white + arab + white_group + arab_group + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata, family="binomial")
out.9controls <- summary(fit.9controls)

fit.9eb <- glm(motive_religious ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial", weights=webal)
out.9eb <- summary(fit.9eb)

## replicate model 10
fit.10 <- glm(motive_political ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")

fit.10controls <- glm(motive_political ~ group + white + arab + white_group + arab_group + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata, family="binomial")
out.10controls <- summary(fit.10controls)

fit.10eb <- glm(motive_political ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial", weights=webal)
out.10eb <- summary(fit.10eb)

## replicate model 11
fit.11 <- glm(motive_insanity ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial")

fit.11controls <- glm(motive_insanity ~ group + white + arab + white_group + arab_group + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata, family="binomial")
out.11controls <- summary(fit.11controls)

fit.11eb <- glm(motive_insanity ~ group + white + arab + white_group + arab_group, data=mydata, family="binomial", weights=webal)
out.11eb <- summary(fit.11eb)

#### Table 5 ####

## replicate model 12
fit.12 <- glm(motive_insanity ~ control + white, data=mydata[which(mydata$all_group==0),], family="binomial")

fit.12controls <- glm(motive_insanity ~ control + white + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata[which(mydata$all_group==0),], family="binomial")
out.12controls <- summary(fit.12controls)

fit.12eb <- glm(motive_insanity ~ control + white, data=mydata[which(mydata$all_group==0),], family="binomial", weights=webal)
out.12eb <- summary(fit.12eb)

## replicate model 13
fit.13 <- glm(motive_insanity ~ group + white_group, data=mydata[which(mydata$all_group==1),], family="binomial")

fit.13controls <- glm(motive_insanity ~ group + white_group + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata[which(mydata$all_group==1),], family="binomial")
out.13controls <- summary(fit.13controls)

fit.13eb <- glm(motive_insanity ~ group + white_group, data=mydata[which(mydata$all_group==1),], family="binomial", weights=webal)
out.13eb <- summary(fit.13eb)


#### June 20 Data ####

#### Table 6 ####

# replicate model 14
fit.14 <- polr(Fgun_policy ~ control + arab, data = mydata2)

fit.14controls <- polr(Fgun_policy ~ control + arab + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata2)

fit.14eb <- polr(Fgun_policy ~ control + arab, weights=webal, data=mydata2)

# model 15
fit.15 <- polr(Fhealth_policy ~ control + arab, data = mydata2)

fit.15controls <- polr(Fhealth_policy ~ control + arab + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata2)

fit.15eb <- polr(Fhealth_policy ~ control + arab, data=mydata2, weights=webal)

# model 16
fit.16 <- polr(Fterror_policy ~ control + arab, data = mydata2)


fit.16controls <- polr(Fterror_policy ~ control + arab + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata2)

fit.16eb <- polr(Fterror_policy ~ control + arab, data=mydata2, weights=webal)

## Table 7 -- most important policies

# Model 17
fit.17 <- glm(gun_mi ~ control + arab, data=mydata2, family="binomial")
summary(fit.17)

fit.17controls <- glm(gun_mi ~ control + arab + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata2, family="binomial")
out.17controls <- summary(fit.17controls)

fit.17eb <- glm(gun_mi ~ control + arab, data=mydata2, family="binomial", weights=webal)
out.17eb <- summary(fit.17eb)

# Model 18
fit.18 <- glm(health_mi ~ control + arab, data=mydata2, family="binomial")
summary(fit.18)

fit.18controls <- glm(health_mi ~ control + arab + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata2, family="binomial")
out.18controls <- summary(fit.18controls)

fit.18eb <- glm(health_mi ~ control + arab, data=mydata2, family="binomial", weights=webal)
out.18eb <- summary(fit.18eb)

# Model 19
fit.19 <- glm(terror_mi ~ control + arab, data=mydata2, family="binomial")
summary(fit.19)

fit.19controls <- glm(terror_mi ~ control + arab + resp_3039 + resp_4049 + resp_5059 + resp_60 + resp_female + resp_eduBach + resp_white + resp_black + resp_hispanic + resp_inc2040 + resp_inc4060 + resp_inc6080 + resp_inc80100 + resp_inc100, data=mydata2, family="binomial")
out.19controls <- summary(fit.19controls)

fit.19eb <- glm(terror_mi ~ control + arab, data=mydata2, family="binomial", weights=webal)
out.19eb <- summary(fit.19eb)



##########
## plot

toPlotData <- function(fit1, fit2, fit3, model) {
    oneD <- function(m, fname) {
        mycoef <- coef(m)
        myci <- confint(m)
        return(data.frame(model=model, Fit=rep(fname,length(mycoef)), Variable=names(mycoef), coef=mycoef, ciL=myci[,1], ciU=myci[,2]))
    }
    
    cleanname <- function(d) {
        d$Variable <- as.character(d$Variable)
        d$Variable[which(d$Variable=="group")] <- "Group Only"
        d$Variable[which(d$Variable=="white")] <- "White"
        d$Variable[which(d$Variable=="white_group")] <- "White Supremacist"
        d$Variable[which(d$Variable=="arab_group")] <- "Islamic Extremist"
        d$Variable[which(d$Variable=="arab")] <- "Arab"
        d$Variable[which(d$Variable=="control")] <- "Control"
        return(d)
    }
    
    temp1 <- oneD(fit1, "Primary")
    temp2 <- oneD(fit2, "Control")
    temp3 <- oneD(fit3, "Entropy Balancing")
    out <- rbind(temp1, temp2, temp3)
    out <- out[which(out$Variable %in% c("control", "group", "white", "arab", "white_group", "arab_group")),]
    out<-cleanname(out)
    return(out)
}

toPlot <- function(tempdata, pd, xlab=TRUE, legend=TRUE) {
    cols <- c("#636363", "#636363", "#636363") # colors
    tempdata$Fit <- factor(tempdata$Fit,levels=c("Primary", "Control", "Entropy Balancing"))
    d <- ggplot(tempdata, aes(x=Variable, y=coef, group=Fit, colour=Fit, linetype=Fit)) + geom_errorbar(aes(ymin=ciL, ymax=ciU), width=.1, position=pd, size=1.1) + scale_linetype_manual(values=c(1,2,3)) + geom_point(position=pd, size=2.5, colour="black") + scale_colour_manual(values=cols)  + labs(x="Variable", y="Effect Size") + theme(axis.title=element_text(size=rel(1.8)), axis.text=element_text(size=rel(1.15)), legend.title=element_text(size=rel(1.8)), legend.text=element_text(size=rel(1.15)), strip.text.x=element_text(size=rel(1.8))) + geom_hline(yintercept=0) + coord_flip()
    if(xlab==FALSE) {
        d <- d + labs(x="")
    }
    if(legend==FALSE) {
        d <- d + theme(legend.position="none")
    }
    return(d)
}

pd1<- toPlotData(fit.1, fit.1controls, fit.1eb, model="Model 1")
pd2<- toPlotData(fit.2, fit.2controls, fit.2eb, model="Model 2")
pd3<- toPlotData(fit.3, fit.3controls, fit.3eb, model="Model 3")
pd4<- toPlotData(fit.4, fit.4controls, fit.4eb, model="Model 4")
pd5<- toPlotData(fit.5, fit.5controls, fit.5eb, model="Model 5")
pd6<- toPlotData(fit.6, fit.6controls, fit.6eb, model="Model 6")
pd7<- toPlotData(fit.7, fit.7controls, fit.7eb, model="Model 7")
pd8<- toPlotData(fit.8, fit.8controls, fit.8eb, model="Model 8")
pd9<- toPlotData(fit.9, fit.9controls, fit.9eb, model="Model 9")
pd10<- toPlotData(fit.10, fit.10controls, fit.10eb, model="Model 10")
pd11<- toPlotData(fit.11, fit.11controls, fit.11eb, model="Model 11")
pd12<- toPlotData(fit.12, fit.12controls, fit.12eb, model="Model 12")
pd13<- toPlotData(fit.13, fit.13controls, fit.13eb, model="Model 13")
pd14<- toPlotData(fit.14, fit.14controls, fit.14eb, model="Model 14")
pd15<- toPlotData(fit.15, fit.15controls, fit.15eb, model="Model 15")
pd16<- toPlotData(fit.16, fit.16controls, fit.16eb, model="Model 16")
pd17<- toPlotData(fit.17, fit.17controls, fit.17eb, model="Model 17")
pd18<- toPlotData(fit.18, fit.18controls, fit.18eb, model="Model 18")
pd19<- toPlotData(fit.19, fit.19controls, fit.19eb, model="Model 19")

# setting a default for dodge
pd <- position_dodge(0.5)

## Figure A1
tempdata <- rbind(pd1, pd2)
d <- toPlot(tempdata, pd)
d <- d + facet_wrap(~ model)
ggsave("appendixA1.pdf", d, width=14, height=7)

## Figure A2
tempdata <- rbind(pd3, pd4)
d <- toPlot(tempdata, pd) + facet_wrap(~ model)
tempdata <- rbind(pd5, pd6)
d2 <- toPlot(tempdata, pd) + facet_wrap(~ model)
ggsave("appendixA2.pdf", arrangeGrob(d, d2, ncol=1), width=14, height=7)

## Figure A3
tempdata <- rbind(pd7, pd8, pd9, pd10, pd11)
tempdata$model <- factor(tempdata$model,levels=c("Model 7", "Model 8", "Model 9", "Model 10", "Model 11"))
d <- toPlot(tempdata, pd)
d <- d + facet_wrap(~ model)
ggsave("appendixA3.pdf", d, width=14, height=7)

## Figure A4
d <- toPlot(pd12, pd, legend=FALSE) + facet_wrap(~ model)
d2 <- toPlot(pd13, pd, xlab=FALSE) + facet_wrap(~ model)
ggsave("appendixA4.pdf", arrangeGrob(d, d2, ncol=2), width=14, height=7)

## figure A5
tempdata <- rbind(pd14, pd15, pd16)
d <- toPlot(tempdata, pd)
d<- d + facet_wrap(~ model)
ggsave("appendixA5.pdf", d, width=14, height=7)

## Figure A6
tempdata <- rbind(pd17, pd18, pd19)
d <- toPlot(tempdata, pd)
d <- d + facet_wrap(~ model)
ggsave("appendixA6.pdf", d, width=14, height=7)

